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ABSTRACT 

Estimation  of  the  phase  of  a  sinusoid  is  an  important  prob¬ 
lem  in  signal  processing.  The  usual  maximum  likelihood 
estimator  is  biased  and  so  can  produce  poor  results,  espe¬ 
cially  at  low  signal-to-noise  ratios  and/or  short  data  records. 
It  is  proven  that  no  unbiased  estimator  exists;  based  on  the 
proof,  a  means  of  obtaining  estimators  with  less  bias  than 
the  maximum  likelihood  estimator  are  proposed. 


the  eigenvalues  and  eigenfunctions  of  the  integral  operator 
and  solve  (1)  via  a  series  solution  [6]. 

2.  PROBLEM  STATEMENT 

We  consider  the  estimation  of  the  phase  0  of  a  sinusoid  em¬ 
bedded  in  noise  or 


1.  INTRODUCTION 

The  phase  of  a  sinusoid  embedded  in  white  Gaussian  noise 
is  a  well  studied  problem.  It  is  usually  solved  by  employ¬ 
ing  the  maximum  likelihood  estimator  (MLE)  [1].  Unfortu¬ 
nately,  the  MEE  is  biased  and  so  can  lead  to  poor  results  at 
low  signal-to-noise  ratios  and/or  short  data  records.  Prob¬ 
lems  such  as  cycle  skipping  in  phase-locked  loops  [2]  are 
directly  attributable  to  this  bias.  Other  practical  problems  of 
interest,  in  addition  to  communications,  are  in  frequency  es¬ 
timation  via  fast  methods  [3]  and  bearing  estimation  in  line 
arrays  [4].  All  of  these  encounter  difficulties  due  to  a  “phase 
wraparound”,  which  is  equivalently  characterized  as  an  es¬ 
timation  bias  error.  In  this  paper  we  investigate  whether 
an  unbiased  phase  estimator  exists  and  if  not,  the  extent  to 
which  a  “nearly  unbiased”  estimator  can  be  implemented. 

A  desirable  property  of  a  statistical  point  estimator  of  a 
parameter  9  is  that  on  average  it  yields  the  true  value  of  the 
parameter.  Specifically,  given  a  random  sample  x  = 
{xo,xi,...  ,a;Ar_i),  where  p(x;0),0  €  0  is  the  probabil¬ 
ity  density  function  (pdf)  of  x,  then  the  statistical  estimator 
(5(x)  of  the  parameter  6  is  said  to  be  unbiased  if  the  expected 
value  of  (5(x)  produces  9',  namely 

ES{tc)  =  j  S{Tc)p{TC]9)dx  =  9,  9  £  Q.  (1) 

A  number  of  well  known  techniques  are  available  to 
solve  integral  equations  in  the  form  of  (1),  such  as  differ¬ 
entiating  or  using  Eourier,  Eaplace  or  Mellin  integral  trans¬ 
forms  [5].  An  alternate  approach  used  here  is  to  determine 


x[n]  =  A  cos[a;on  4-  0]  4-  w[n],  n  =  0, 1, . . .  ,  N  —  1 

where  w[n]  is  white  Gaussian  noise  with  unknown  variance 
(T^ ;  furthermore,  the  amplitude  A  is  considered  unknown 
while  the  frequency  uJq  is  assumed  to  be  known.  All  the 
information  regarding  the  parameter  9  is  summarized  by  the 
jointly  sufficient  statistics  given  by  [1]: 


N-l 

Ti  =  x[n]  cos[a;on], 

n— 0 
N-l 

T2  =  x[n]  sin[a;on]. 

n=0 


Consequently,  all  inference  can  be  based  on  the  joint 
pdf  of  T  =  (Ti ,  T2)  and  is  given  by: 


where 


h{9)  =  exp 


N-l 


cos^[a;on  4-  9] 


n=0 
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and 


E  = 


N-l 

cos^[na;o] 

n=0 

N-l 

cos[na;o]  sin[na;o] 

n=0 


N-l 

cos[na;o]  sin[na;o] 

n=0 

N-l 

sin^[na;o] 

n=0 


is  a  27r  periodic,  symmetric  kernel.  Furthermore  using  [8, 
3.462.5]  the  kernel  can  be  expressed  in  terms  of  $(a;),  the 
standard  error  function  [8,  8.250.1]: 


K{e  -  ct>) 


m 

TT 


is  the  covariance  matrix. 


where 


We  seek  an  unbiased  estimator  (5(T)  of  9\  namely 


h  = 


aVn 

a 


cos(0  —  (f)). 


/oo  roo 

/  S(T)p(T;0)dT  =0,  -tt  <  6»  <  tt.  (2) 

-OO  J  — oo 

For  mathematical  simplicity,  it  is  assumed  that  lOq  =  ; 

consequently,  E  =  yl,  |  det  E|^  =  -y 
and  therefore 


3.  INTEGRAL  EQUATIONS 

Integral  equations  are  equations  in  which  an  unknown  func¬ 
tion  appears  under  the  integral  sign.  The  relevant  form  con¬ 
sidered  here  is  the  Fredholm  equation  of  the  first  kind  [7] 
with  a  27r  periodic,  symmetric  kernel  K{x,  t)  \  namely. 


P(T;  0)  =  exp  ( ^  (Ti  cos[6»]  -  T2  sin[6»])^ 

TTa^JN  \a^  ) 

X  exp(^-^(T2+T|)j  ,  (3) 

with 

M^)=exp(-^)- 

Transforming  (3)  to  polar  coordinates  via  Ti  =  pcos[^]  and 
T2  =  —p  sin[^] ,  where  the  latter  minus  sign  has  been  intro¬ 
duced  in  order  to  subsequently  obtain  a  symmetric  kernel, 
yields 


f{x)=  K{x  -  t)g{t)dt,  K{x,t)  =  K{t,x)  (6) 

J  —TT 

where  f{x)  is  given  and  g{x)  is  an  unknown  function.  A 
basic  result  [6]  is  that  given  a  continuous  real  and  symmet¬ 
ric  kernel  and  a  continuous  f{x),  then  a  solution  to  (6)  ex¬ 
ists  only  if  the  given  function  f{x)  can  be  expanded  in  se¬ 
ries  whose  basis  functions  are  orthonormal  eigenfunctions, 
'^k{x)  of  the  kernel  K{x  —  t)  \  namely, 

OO 

f{x)=  Y.  fk^k{x)  (7) 

k=  —  oo 

with  coefficients 

fk=[  f(t)'^k(t)dt. 

J  —TT 

Next,  expanding  the  unknown  function  g{x)  in  a  Fourier 
series  in  terms  of  ^i.{x)  \  namely. 


Furthermore,  it  is  assumed  that  the  estimator  is  scale  invari¬ 
ant;  namely,  the  same  estimate  is  obtained  if  T  is  multiplied 
by  any  constant  c  >  0.  Consequently,  S{p,(f>)  will  depend 
only  on  f.  Under  this  invariance  assumption,  the  polar  form 
of  (2)  is  given  by  the  following  integral  equation: 


OO 

9{x)  =  Y  Sk'^kix) 

k=  —  oo 

and  substituting  into  (6)  we  have  that 


S{(t>)K{0  -  <f)df  =  0, 


— TT  <  9  <  TT 


(5) 


K{x  —  f)g{f)dt 


where 


K{0  -  f) 


h{0) 


pexp 


cos{9  —  (f) 


dp 


/TT  ^ 

K{x-t)  Y  9k'^k(t)dt 

-TT  1.  -- 


k=  —  oo 


Y  9k  _ 

k=—oo 


K{x 


f)'9k{t)dt 


OO 

Y  9k^k'^k{x)  (8) 

k=  —  oo 


where 


A*  =  /  K{x  -t)'^k{'t)dt. 

J  —TT 

Examples  of  the  kernel  K{t)  for  values  of  /3  =  are 
shown  in  Figure  1 .  Equating  (7)  with  (8),  it  can  be  seen  that 
the  solution  to  (6)  can  be  expressed  as  a  series  in  terms  of 
the  Fourier  series  coefficients  of  f{x)  and  the  eigenvalues 
and  eigenfunctions  of  the  kernel  [6] ;  namely 


Furthermore,  the  complex  Fourier  series  representation  of  9 
is 

OO 

6  =  ^  6k  exp(^A:0) 

k=  —  oo 

where  the  complex  Fourier  coefficients  6k  satisfy  the  rela¬ 
tionship 


provided 


9{x)  = 


E 


/* 

A* 


'^k(x) 


E 


A* 


<  OO. 


(9) 


6k 


1  r 
^  J- 


t  exp(—tkt)dt 


»(-!)*■ 

k 

0 


k  =  Q. 


Therefore  from  (9),  the  series  expression  for  the  unbiased 
estimate  of  phase  is  given  by 


1.4| 


-  -  p  =  0  dB 
-  p  =  6dB 
p=12dB 


d{4>)=  E  ^exp(^A:^) 


(10) 


0.6 1 


0.4 1 


provided  I  I  <00. 

The  Fourier  coefficients  9k  explicitly  decay  as  Ijk.  In 
general,  the  rate  of  decay  of  Fourier  coefficients  is  related  to 
the  smoothness  of  the  function  [9].  In  particular,  the  deriva¬ 
tives  of  the  kernel  with  respect  to  6  exist  for  all  orders  and 
it  is  easily  shown  that 


|Ai;|  - 


< 


1 

1 


d^K{t) 

dt^ 

d^K{t) 

dt^ 


exp(tkt)dt 
exp(tkt)dt . 


Fig.  1.  Integral  equation  kernel. 

We  note  that  the  kernel  is  an  analytic  and  periodic  function 
of  6',  consequently,  this  kernel  can  generate  only  analytic 
and  periodic  functions  of  6,  regardless  of  S{(f>).  However, 
the  desired  right  hand  side  (RHS)  of  (5);  namely  6,  is  an¬ 
alytic  but  not  periodic.  Conversely,  if  RHS  is  periodically 
extended,  that  extension  is  not  analytic.  Consequently,  no 
exact  solution  exists.  This  is  rigorously  proven  by  demon¬ 
strating  that  a  series  solution  does  not  exist. 

It  is  easily  shown  [7]  that  the  eigenfunctions  of  K  {6  —  (f>) 
are  complex  exponentials,  exp{ik(f>)  and  consequently  the 
eigenvalues  are 


This  implies  that  Xk  decays  faster  than  1/A:;  consequently 

OO  I  i?  |2 


and  therefore  no  unbiased  estimate  of  phase  S{(f>)  exists. 

4.  MAXIMUM  LIKELIHOOD  ESTIMATION 

The  maximum  likelihood  estimate  of  sinusoidal  phase  is 
given  by  [  1  ] 


A*  —  /'  K(t)  exp(—tkt)dt . 


6  =  —  arctan 


(11) 


Applying  the  expectation  operator  we  obtain 


Ee 


ep(Ti,T2-,e)dTidT2 
—  arctan  p{Ti,T‘2',6)dTidT‘2  . 


Transforming  this  expression  to  polar  coordinates  via  Ti  = 
pcos[^]  andT2  =  — psin(^)  yields 


E§ 


4>P{P,  4>',  ^)dp  dcj) 


—  TT 


Jo 

<i>K{6  —  <i>)d<i> . 


(12) 


Expanding  (p  in  terms  of  the  eigenfunctions  of  K (9  —  (f))\ 
namely,  a  Fourier  series  and  substituting  into  (12): 


s 

I  0 


-  - p  =  0dB 

-  p=6dB 
P  =  12dB 


Fig.  2.  Expected  maximum  likelihood  estimate. 


E§ 


/TT  ^ 

4>k  exp(tk4>)K(9  —  4>)d4> 

k—  —  oo 

OO 

4>k  /  exp(tk4>)K(9  —  4>)d4> 

_  J  —If 


k—  —  oo 


OO 

X!  ('i;A^;exp(^A:6»)  (13) 

k=  —  oo 


where  we  have  used  the  fact  that  ([>1.  =  9k- 

It  is  seen  that  the  expected  value  of  the  MEE  depends  on 
the  reciprocal  of  the  eigenvalues  of  S{(f>)  and  is  the  source  of 
the  MEE’s  bias.  Examples  of  the  average  MEE,  illustrating 
this  bias,  are  shown  in  Figure  2.  Future  work  will  consider 
estimators  of  the  form 


k—  —  oo  ^ 

where  the  weights  Wk  are  the  Fourier  coefficients  of  a  win¬ 
dow  function  w.  The  weights  are  chosen  such  that  the  es¬ 
timator  obtained  from  the  convolution  of  w  with  9  has  less 
bias  than  the  MEE,  at  the  possible  expense  of  variance. 


5.  CONCLUSIONS 
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